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ABSTRACT 

An  iterative  approximation  algorithm  is  presented  for  calculating  the  stationary  queue  size 
probabilities  of  tandem  queueing  systems  subject  to  blocking.  The  algorithm  combines  a  decom¬ 
position/aggregation  technique  with  exact  analytical  results  for  two  node  systems.  It  applies  to 
tandem  lines  where  failure  type  servers  with  phase-type  service  and  repair  time  distributions  are  in 
attendance  under  various  blocking  disciplines 
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1.  INTRODUCTION 


Consider  a  production  system  composed  of  many  processing  stages  (servers)  through  which  the 
parts  (jobs)  must  pass  in  some  prespecified  order.  In  practice,  the  servers’  behavior  is  modulated 
by  several  factors,  such  as  fluctuations  in  the  service  times  and  interruptions  of  service  due  to  server 
breakdowns.  Although  this  service  variability  greatly  affects  the  performance  of  such  a  system,  its 
influence  can  be  somewhat  mitigated  by  using  intermediate  storage  spaces  (buffers)  between  the 
servers.  However,  because  of  physical  limitations  on  the  buffer  sizes,  the  flow  of  jobs  through  the 
system  may  get  blocked. 

The  blocking  systems  considered  here  are  composed  of  a  series  configuration  of  service  stations 
with  finite  intermediate  buffers  between  the  stations.  A  typical  tandem  configuration  is  shown  in 
Figure  1.  At  each  station  some  work  is  performed  on  a  job  which  is  then  passed  on  to  the  next 
station;  upon  completion  of  service  at  the  last  station,  the  job  immediately  leaves  the  system. 
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Figure  1. 

The  literature  considers  two  distinct  blocking  policies  for  transfer  lines  with  finite  buffers.  In 
order  to  describe  them,  let  S  be  one  of  the  servers  in  Figure  1.  Under  the  first  policy,  called 
immediate  blocking ,  at  a  time  of  service  completion,  the  server  S  is  blocked  if  the  downstream 
buffer  becomes  full  as  a  result  of  this  service  completion;  the  server  S  remains  blocked  until  a  space 
becomes  available  in  the  downstream  buffer,  at  which  time  the  server  resumes  service  and  begins 
processing  its  next  job  (if  any).  Under  the  second  policy,  called  non-immediate  blocking,  the  server 
S  is  blocked  at  a  service  completion  time  if  the  job  that  has  just  completed  service  finds  the  next 
buffer  to  be  full.  As  soon  as  a  space  becomes  available  in  the  next  buffer  this  job  moves  downstream 
without  receiving  any  further  service  at  server  S,  which  then  begins  processing  its  next  job  (if  any). 
As  noted  by  Altiok  and  Stidham  [2],  these  two  blocking  policies  are  in  general  not  equivalent  in 
that  the  solution  of  the  system  under  one  policy  cannot  be  obtained  from  the  solution  of  the  system 
(with  possibly  different  parameter  values)  under  the  other  policy;  such  an  equivalence  does  however 
hold  for  two  station  systems. 

Queueing  networks  with  blocking  have  been  studied  by  researchers  from  different  research 
communities,  as  evidenced  by  the  annotated  bibliography  on  blocking  systems  compiled  in  [7].  A 
basic  assumption  made  in  all  the  models  studied  in  the  surveyed  literature  is  that  the  last  stage  is 
never  blocked,  i.  e.,  there  is  always  space  available  for  a  job  whose  service  has  been  completed  at 
the  last  server. 

Queueing  networks  with  blocking  are  typically  very  difficult  to  analyze,  and  closed  form  results 
at  steady  state  (or  otherwise)  are  usually  not  available.  Although  a  rigorous  Markovian  analysis  was 
performed  early  on  by  Hunt  [9],  to  date  strikingly  few  results  have  been  obtained,  with  the  bulk  of 
the  work  focusing  on  continuous-time  models.  This  state  of  affairs  points  to  the  need  of  developing 
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efficient  approximation  techniques  to  evaluate  performance  measures  of  interest.  In  this  paper,  an 
iterative  approximation  scheme  is  presented  for  finding  the  steady-state  marginal  probabilities  of 
the  queue  sizes  in  tandem  queueing  systems  with  finite  capacity  buffers  and  phase-type  (PH-type) 
servers.  The  algorithm  is  based  on  the  analytical  results  obtained  in  [8]  for  two  node  tandem 
systems,  and  is  presented  under  the  immediate  blocking  policy.  The  same  approximation  scheme 
also  applies  to  tandem  systems  under  the  non-immediate  blocking  policy  in  view  of  the  above- 
mentioned  equivalence  for  two  node  systems,  and  is  only  briefly  outlined  in  the  interest  of  brevity. 
In  light  of  the  results  derived  in  [8],  the  proposed  approximation  scheme  is  also  applicable  to  systems 
with  failure  type  servers  with  PH-type  service  and  repair  distributions,  and  to  systems  where  jobs 
serviced  at  the  ( i  +  l)st  server  can  be  fed  back  to  the  ith  buffer.  Since  the  approximation  scheme  is 
based  on  the  effective  solution  of  two  node  systems,  the  relevant  results  of  [8]  are  summarized  in  the 
Appendix.  The  reader  is  also  referred  to  the  work  of  Buzacott  and  Kostelski  [4]  for  an  algorithmic 
solution  of  two  node  production  systems  with  two-stage  Coxian  service  distributions. 

Most  approximation  algorithms  use  the  flow  conservation  principle  to  decompose  the  tandem 
model  into  simple  two  node  models.  The  approximation  reported  here  uses  this  viewpoint  in  the 
form  taken  by  Altiok  in  [1],  and  represents  the  effective  service  distribution  of  a  server  by  considering 
all  the  servers  upstream  of  this  server  in  order  to  capture  the  effect  of  blocking.  However,  the 
algorithm  presented  here  differs  from  Altiok’s  in  that  instead  of  approximating  the  arrivals  to  the 
ith  buffer  by  a  Poisson  process,  another  similar  effective  representation  is  introduced  to  capture 
the  effect  of  idling.  Once  this  decomposition  step  is  taken,  approximations  are  made  to  express  the 
effective  representations  iteratively  in  terms  of  quantities  that  can  be  obtained  by  solving  two  node 
models. 

Although  the  discussion  is  given  for  the  discrete-time  formulation,  the  ideas  presented  here 
apply  mutatis  mutandis  to  the  continuous-time  formulation.  The  modeling  of  time  as  a  discrete 
parameter  could  be  motivated  by  the  fact  that  service  times  in  manufacturing  systems  are  often 
constant,  the  main  source  of  randomness  being  introduced  by  the  possibility  of  server  failures. 
The  accuracy  of  the  algorithm  is  validated  through  numerical  examples,  both  for  continuous  and 
discrete-time  systems.  Comparison  of  the  results  against  simulations  indicates  reasonable  accuracy 
under  both  blocking  policies  even  in  the  presence  of  significant  blocking.  The  algorithm  is  also 
compared  to  another  approximation  algorithm  proposed  by  Jun  and  Perros  [10]  for  open  tandem 
queueing  systems  with  two-stage  Coxian  service  distributions. 

Phase  Type  Distributions 

To  fix  the  notation  and  terminology,  it  is  convenient  at  this  point  to  briefly  describe  the  class 
of  discrete-time  PH-type  distributions.  The  reader  is  invited  to  consult  the  monograph  by  Neuts 
[11]  for  additional  information  on  this  topic  and  for  a  similar  development  on  continuous-time  PH- 
distributions.  The  following  notation  is  used  hereafter:  For  any  positive  integer  r,  the  r  x  r  identity 
matrix  is  denoted  by  Ir  and  the  r  X  1  column  vector  of  ones  is  denoted  by  er,  while  the  1  X  r 
dimensional  row  vector  with  all  zero  entries  is  denoted  by  0r. 

Loosely  speaking,  a  discrete-time  PH-type  distribution  is  any  probability  distribution  on  the 
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non-negative  integers  which  can  be  realized  as  the  distribution  of  the  time  until  absorption  in  a 
discrete-time  finite  state  space  Markov  chain  with  a  single  absorbing  state.  The  class  of  PH-type 
distributions  includes  well-known  distributions  such  as  the  Generalized  Negative  Binomial  and 
Hypergeometric  distributions,  as  well  as  any  distribution  with  finite  support  on  the  non- negative 
integers. 

Specifically,  consider  a  discrete-time  Markov  chain  on  the  state  space  {1,2,. .  .,m+  1),  where 
the  states  {l,...,m}  are  transient  and  the  state  m  +  1  is  absorbing.  The  chain  starts  in  state  j 
with  probability  aj  and  evolves  according  to  the  matrix  P  of  one-step  transition  probabilities  with 


P  = 


(a,am+i)  and  a  =  (ax,a2,. .  .,am)  , 


(1.1) 


where 

Q  is  a  m  X  m  substochastic  matrix, 

p  is  a  m  x  1  column  vector  of  absorption  probabilities  into  the  absorbing  state  m  +  1  from 
the  transient  states  {1, . .  .,m}, 
a  is  a  1  X  m  row  vector  of  initialization  probabilities. 

The  PH-type  probability  distribution  associated  with  the  pair  (a,  Q)  is  the  distribution  function 
F  on  the  non-negative  integers  with  probability  mass  function  {qk,  k  =  0, 1, . . .}  given  by 

(  &m- 1-1  j  k  ~  0  , 

qk  =  (1.2) 

l  aQk  lp  ,  k  >  1  . 


The  pair  (a,<2)  is  called  the  representation  of  the  distribution  F,  which  it  uniquely  determines. 
The  converse  is  not  true  in  that  a  given  PH-type  distribution  function  F  admits  infinitely  many 
PH-representations. 

The  following  probabilistic  construction  is  given  in  Neuts  [11,  p.  48].  Upon  absorption  into 
state  m  +  1,  independent  multinomial  trials  with  probabilities  (a,am+i)  are  instantaneously  and 
repeatedly  performed  until  one  of  the  alternatives  1,2,  ...,m  occurs,  say  j.  The  process  is  now 
restarted  in  the  state  j  and  the  same  procedure  is  repeated  at  the  next  absorption.  Upon  indefinitely 
continuing  this  procedure,  a  new  Markov  process  is  constructed  on  {1, 2, ... ,  m}  with  the  state  m+ 1 
as  an  instantaneous  state,  and  with  matrix  Q*  of  one-step  transition  probabilities  given  by 


Q*  =  Q  +  r - - - V  a  • 

1  Q’rn+l 


(1.3) 


The  representation  (a,  Q )  is  said  irreducible  if  the  Markov  chain  on  {1, . . . ,  m}  with  one-step  tran¬ 
sition  matrix  Q*  is  irreducible.  Every  PH-type  distribution  function  F  admits  an  irreducible  rep¬ 
resentation  [11,  p.  49]. 
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2.  THE  DISCRETE-TIME  MODEL 


Consider  the  tandem  system  shown  in  Figure  1,  which  consists  of  N  nodes  with  N  —  1  interme¬ 
diate  finite  capacity  buffers  of  sizes  A',,  1  <  i  <  N,  including  the  jobs  being  served  by  the  servers. 
Each  node  is  attended  by  a  single  server  which  operates  according  to  the  FCFS  (first-come-first- 
served)  queueing  discipline.  For  1  <  i  <  N,  the  service  times  at  the  ith  node  are  assumed  to  be 
independent  and  identically  distributed  (i.i.d)  with  common  discrete  PH-distribution  Fi  given  by 
the  irreducible  PH-representation  («;,  Qi)  of  order  m;.  The  m;  X  1  column  vectors  of  absorption 
(service  completion  at  the  ith  node)  probabilities  is  denoted  by  pi ,  and  in  order  to  avoid  cases  of 
limited  interest,  the  vector  of  initialization  probabilities  is  assumed  to  satisfy  a;emj  =  1.  The 
matrix  (Imj  —  Qi)  is  assumed  nonsingular ,  so  that  the  service  completion  from  any  initial  phase  is 
certain  [11,  p.  45].  The  service  times  at  different  servers  are  also  assumed  mutually  independent. 
It  is  assumed  that  the  last  server  is  never  blocked  while  the  first  server  is  always  busy ,  i.  e.,  there 
is  an  infinite  supply  of  exogenous  jobs.  The  situation  where  the  first  server  is  generating  arrivals  to 
the  first  buffer  according  to  a  PH-renewal  process  can  also  be  approximated  by  a  similar  analysis. 
Let  the  set  S x  of  service  phases  for  the  ith  server  be  defined  by 

S'  :=  {sj  :  1  <  /  <  m;}  ,  1  <  i  <  N, 

and  define  the  scalars 

i 

1  <i<j<N. 

k=i 

3.  DECOMPOSITION  AND  APPROXIMATION  METHOD 

The  steady-state  marginal  queue  length  distribution  for  the  ith  node  could  be  calculated  from 
Theorem  A.l  of  the  Appendix  if  the  ( i  +  l)st  server  were  not  subject  to  blocking  and  the  ith  server 
were  always  busy.  However,  at  a  service  completion  epoch  for  the  ( i  +  l)s<  server,  this  server  may 
be  blocked  and  may  remain  blocked  until  there  is  a  departure  from  the  ( i  +  l)s<  buffer.  Similarly, 
at  a  service  completion  at  the  ith  server,  this  server  may  become  idle.  The  decomposition  method 
proposed  here  amounts  to  finding  equivalent  PH-representations  for  the  ith  and  the  (i  +  l)s*  servers 
that  incorporate  the  effects  of  idling  and  blocking,  respectively. 

For  l<i,j<N,l<l<mj  and  t  >  0,  define  the  events 

Bl{t)  :  the  ith  server  is  blocked  at  time  t, 

B\'3(t)  :  the  servers  k,  i  <  k  <  j,  are  all  blocked  at  time  t,  and  the  jth  server  is  not  blocked  and 

in  service  phase  sfi 

P(t)  :  the  ith  server  is  idle  at  time  t, 

I\'\t, )  :  the  servers  k,  j  <  k  <  i,  are  all  idle  at  time  t,  and  the  jth  server  is  not  idle  and  in 

service  phase  sj. 
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Assuming  that  the  limits  exist,  define  the  probabilities 


Ph  :=  limn  iWT*)]  , 
P'di  ■=  limn  P[Bi'j(Ti)}  , 
P\  :=  limn  P[P(Ti)}  , 
Pjf  :=  lim„  P{lj’j(Ti)}  , 


1  <  i  <  N, 

1  5=  *  <  j  <  A,  1  <  l  <  mj, 
1  <  i  <  A, 

1  <  <  j  <  A,  1  <  l  <  rrtj, 


(3.1) 


where  { T n  =  1, 2, . . .}  is  the  sequence  of  service  completion  times  at  node  i,  and  lim„  is  taken 
with  n  going  to  infinity. 


Decomposition 

Fix  1  <  i  <  A  and  consider  the  nth  service  completion  time  T^+1  in  the  (z  +  l)st  server. 
This  server  is  either  blocked  with  probability  P[B1+1  (Tf+1)}  or  with  complementary  probability 
1  —  P[Bt+1  (T^+1 )]  it  starts  serving  its  next  job  (if  any)  according  to  the  initialization  vector  ai+i .  In 
the  case  of  blocking,  the  event  il;!+1,J(T^+1)  takes  place  for  some  j,  i+l  <  j  <  N ,  and  /,  1  <  l  <  mj, 
in  which  case  the  effective  service  time  of  the  blocked  job  at  the  (i  -f  l)st  server  is  equal  to  the 
residual  service  time  in  the  jth  server  plus  the  sum  of  the  service  times  of  all  the  blocked  servers 
j  —  1,. .  ,,i  +  2  (time  to  unblock),  plus  the  service  time  in  the  (j  +  l)st  server.  Therefore,  if  the 
( i  +  l)sf  server  is  blocked  at  time  T^+1,  then  the  service  phase  of  the  ( i  +  l)st  server  is  effectively 
reinitialized  at  phase  sj  with  probability  P[j3i!+1’J(T^+1)]  and  the  transitions  occur  within  the  set 
according  to  the  matrix  Qj.  Upon  service  completion  in  the  jth  server,  service  in  the  ( j  —  l)s< 
server  is  initialized  according  to  oij-i,  thus  causing  a  transition  from  the  set  S *  to  the  set  S'J~1  in 
the  effective  representation  of  the  (i  +  l)st  server.  This  argument  can  be  repeated  upon  moving 
upstream  from  station  j  to  station  ( i  +  1)  and  finally,  upon  service  completion  in  the  (z  +  2)nd 
server,  the  blocking  period  of  the  ( i  +  l)st  server  ends  and  this  server  resumes  service  according  to 
the  initialization  vector 

Therefore  in  steady  state,  the  effect  of  blocking  on  the  ( i  +  l)st  server  can  be  captured  by 
considering  a  new  PH-distribution  for  this  server  with  set  of  service  phases  given  by  (Jjlt+i  . 
This  effective  service  time  distribution  of  order  is  obtained  by  considering  all  the  servers 

downstream  from  the  ( i  +  l)s<  server,  and  is  denoted  by  the  pair  (a2(i  +  1),Q2(*  +  1))-  If  the 
corresponding  SfJ.1  X  1  column  vector  of  absorption  probabilities  is  denoted  by  p2(i  +  1),  then 


Q2{i  +  1)  — 

/  Qi+l 

Pi+2  OCi+l  Qi  +  2 

\ 

,  P2(*+l)  = 

/  Pi+ 1  \ 

\ 

PN&N- 1  Qn > 

VOL  / 

and  the  1  X  initialization  probability  vector  a2(i  +  1)  of  this  effective  representation  of  the 

(i  +  l)s<  server  is  given  by 

a2(i  +  l)=  ((l-P)+1)a1+1,P^1’i+2,...,Pi+1’7V) 


(3.26) 


for  1  <  i  <  N .  For  1  <  i  <  j  <  N,  the  1  X  mj  row  vector  has  entries  1  <  /  <  m,j. 

With  this  effective  PH-representation,  the  ( i  +  l)st  server  is  never  blocked,  but  a  service 
completion  is  allowed  only  from  the  phases  that  correspond  to  a  service  in  the  (i  +  l)s<  server.  The 
subscript  2  indicates  that  the  representation  (cM*  +  1)>  CM*  +  1))  is  the  effective  representation  for 
the  second  node  server  when  considering  the  ith  buffer. 

A  similar  argument  can  now  be  developed  to  capture  the  effect  of  idling,  in  which  case  for 
1  <  i  <  A ,  the  effective  service  of  the  ith  server  can  be  represented  by  a  PH-distribution  of 
dimension  E)  with  set  of  service  phases  (J‘=1  5J.  The  effective  representation  of  the  ith  server  as 
the  first  node  server  in  the  two  node  equivalent  representation  when  considering  the  ith  buffer  is 
denoted  by  (oi(i),  CM*))-  It  is  given  by 


Qi(0  = 


(  Qi 

Pi-i  a;  Qi-i 


V 


\ 


(  Pi  \ 
1  0T  ' 


>  Pi(*)  = 


1 


Pi  <*i  Q\  / 

where  Pi(f)  is  the  X  1  column  vector  of  absorption  probabilities,  and 

«!(*•)=  ((l-piKpr-1,...^;-1) 


Vo f  / 


(3.3  a) 


(3.3 b) 


for  1  <  i  <  N.  The  1  x  mj  row  vector  P\’*  has  entries  P\ 1  <  l  <  mj  and  1  <  j  <  i  <  N. 

Therefore,  by  considering  the  effective  representations  (3.2)  and  (3.3)  for  each  buffer  i,  1  < 
i  <  iV,  the  tandem  system  in  Figure  1  can  be  decomposed  into  N  —  1  two  node  systems  of  the  form 
shown  in  Figure  2. 


MO.  Qitt))^0!  KJ  l®->2(i  +  1),  Qfii  +  1)) 


Figure  2. 

Since  the  matrices  Imi  —  Qi,  1  <  i  <  N ,  are  all  invertible,  so  are  the  matrices  —  Q\{i)  and 
—  Q 2(1  +  1),  1  <  i  <  N .  Furthermore,  since  the  first  (resp.  last)  node  server  of  Figure  1  is 
never  starved  (resp.  blocked),  the  first  (resp.  second)  node  server  of  Figure  2  is  also  never  starved 
(resp.  blocked),  and  the  solution  technique  given  in  [8]  thus  applies. 

Approximations 

Although  the  tandem  line  of  Figure  1  has  been  decomposed  into  N  —  1  two  node  systems 
of  the  type  displayed  in  Figure  2,  these  systems  are  coupled  through  the  initialization  vectors 
k  =  1,2  and  1  <  i  <  N,  and  approximations  are  therefore  needed  in  order  to  compute  the 
vectors  a^( i). 

As  in  [8],  a  Markov  chain  is  associated  with  each  one  of  the  two  node  models  shown  in  Figure 
2.  For  1  <  i  <  N,  let  the  vector  u(f)  (resp.  v(i))  be  the  steady-state  probability  vector  of  this 
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Markov  chain  embedded  at  arrival  (resp.  departure)  epochs  at  the  first  (resp.  second)  node.  The 
1  X  Ei  block  component  «o(0  of  v(i )  (see  Appendix)  is  the  steady-state  probability  vector  that  the 
second  server  becomes  idle  after  a  service  completion.  If  k  =  E)  +  l  ,  1  <  j  <  i,  1  <  l  <  rrij-i, 
then  the  kth  component  of  va(i)  is  the  joint  probability  that  at  a  service  completion  epoch  the 
second  server  is  idle  and  the  phase  of  the  first  node  server  is  sj-1.  In  view  of  (3.1)  and  of  the 
construction  of  Qi(i),  the  kth  component  of  the  vector  wo(*)  approximates  the  Ith  component  of 
the  vector  Pj+1,J_1,  i.  e.,  the  steady-state  probability  that  at  a  service  completion  epoch  of  the 
(i  -f  l)s*  server,  the  buffers  i,i  —  1, . . .,  j  —  1  are  all  empty  and  the  (j  -  l)th  server  is  in  s j~*  service 
phase. 

Similarly,  if  k  =  Eg_i  +  l  for  i  <  j  <  N  and  1  <  l  <  1,  then  the  kth  component  of  the 

1  x  Eg.!  block  «/f(-i(t)  of  u(i)  is  the  joint  steady-state  probability  that  at  a  service  completion 
epoch  in  the  first  server,  it  is  blocked  and  the  phase  of  the  second  node  server  is  Sj+1.  The  kth 
component  of  the  vector  ttK{-i(*)  therefore  approximates  the  Ith  component  of  the  vector  PgJ+1, 
i.  e.,  the  steady-state  probability  that  at  a  service  completion  epoch  of  the  ith  server,  the  servers 
i,i  +  1, . .  .,j  are  all  blocked  and  the  (j  -f  l)s<  server  is  in  its  Sj+1  service  phase. 

In  vector  form  these  approximations  can  be  rewritten  as 

. . . ,  p;+1-x)  «  vo(i)  ,  1  <  i  <  N  (3.4a) 

and 

(P^+1, . . . ,  P^N)  «  u(/c,._1)(0  ,  1  <  *  <  N.  (3.46) 

The  scalars  Po(«)  and  Pp( i)  defined  by  the  relations 

Po(0  =  uo(*)  es<i  and  PF(i )  =  «(/<-._!)(*)  ,  1  <i  <  N, 

thus  provide  approximations  to  Pj+1  and  PXB ,  respectively.  Consequently,  the  initialization  vectors 
«! (i)  and  ct2 (i  +  1)  in  (3.3b)  and  (3.2b)  can  be  approximated  by  the  quantities  a“ (i)  and  a%(i  +  1) 
given  by 

tti(0  =  [  C1  -  Po(i  ~  1)  )  ai,v0(i  -  1)]  (3.5a) 

and 

a%(i+  1)  =  [(1  -  PF(i  +  1))  ai+i  ,w/ci+1_i(i  +  1)]  .  (3.56) 

This  approximation  leads  to  an  iterative  approach  based  on  the  two  node  analysis  summarized 
in  the  Appendix.  The  marginal  probability  distribution  of  the  queue  sizes  at  buffer  i  can  be  obtained 
once  the  vectors  vo(i  —  1)  and  U(/,:.+1_i)(i  +  1)  are  known.  These  vectors  are  calculated  from  the 
two  node  approximations  of  the  ( i  -  l)si  and  ( i  +  l)sf  queues,  respectively.  More  precisely,  with 
some  abuse  in  the  notation,  denoting  the  iteration  count  by  a  superscript,  the  steps  of  the  iterative 
procedure  can  be  described  as  follows: 
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Step  1.  Select  an  initial  value  for  the  vectors  _!)(*')  for  all  1  <  i  <  N  and  compute  a°2(i)  from 
(3.5b).  Set  a2(N)  =  a jy. 

Step  2.  At  iteration  n; 

For  i  =  1,2 ..,N  —  1, 

Starting  from  the  first  queue,  solve  the  effective  two  node  systems  for  each  buffer: 

a.  Set  (1)  =  «i  and  a2(N)  =  aN. 

b.  Compute  «?(*)  and  a?(i)  by  using  v£(i  -  1)  and  u^_r)(i),  respectively. 

c.  Obtain  the  vectors  v£(i)  and  «“K._i)(*)  bY  using  the  results  of  the  Appendix. 

Step  3.  Test  for  a  convergence  criterion,  namely  that  all  the  marginal  queue  length  probabilities 
in  successive  iterations  are  required  to  be  within  e  of  each  other.  If  this  requirement  is 
not  satisfied,  set  n  to  n  -f  1  and  go  to  Step  2. 

Although  no  proof  of  convergence  is  available,  the  algorithm  has  been  tested  on  many  examples 
and  has  always  converged  with  e  =  10-4  in  about  only  5  iterations.  Since  at  each  node  all  the 
servers  are  taken  into  account,  the  computational  complexity  of  the  algorithm  grows  quadratically 
with  N,  the  number  of  servers  in  tandem.  In  fact,  the  proposed  algorithm  can  be  computationally 
prohibitive  for  very  long  tandem  lines  with  high  dimensional  PH-type  distributions.  However,  in 
such  cases,  at  the  expense  of  introducing  further  approximations,  the  following  modifications  can 
easily  be  incorporated.  Firstly,  the  computational  complexity  of  the  algorithm  can  be  made  linear  in 
N  by  considering  only  a  few  immediate  neighboring  nodes  in  the  effective  representations.  Secondly, 
the  effective  representations  Qk(i),  k  =  1,2  and  1  <  i  <  N,  can  be  approximated  by  lower  order 
PH-distributions.  On  the  other  hand,  the  algorithm  is  eminently  suitable  for  parallel  computations 
and  can  be  made  much  faster  by  implementing  it  on  a  parallel  machine.  With  a  different  processor 
assigned  to  the  solution  of  each  two  node  system,  processor  i  needs  information  only  from  its 
immediate  neighbors,  namely  from  processors  i  +  1  and  i  —  1,  to  update  the  initialization  vectors 
cci(i)  and  ao(i  +  1)- 

4.  DECOMPOSITION  FOR  THE  NON-IMMEDIATE  BLOCKING  POLICY 

The  decomposition  described  in  Section  3  for  systems  that  operate  under  the  immediate  block¬ 
ing  strategy  can  also  be  made  for  systems  operating  under  the  non-immediate  blocking  strategy. 
Since  only  the  type  of  blocking  is  different,  the  representation  (3.3)  given  for  the  pair  («i(z),Qi(0) 
remains  unchanged. 

However,  by  an  argument  similar  to  the  one  given  in  Section  3,  the  effective  representation  for 
the  pair  (a2(i  +  1),  Q2(i  +  1))  is  now  seen  to  be 


Q?(i  +  1)  — 


(  Qi  + 1  Pi+ 1  (P j 
Qi+ 2 


i+l,;+2 


Qn- 


Pb+1'N)\ 


Qn  ) 


,  p2(i+  1)  = 


/(l-PiB+l)Pi+l\ 

Vi+2 


PN 


(4-1 a) 
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and 

£*2(2  d  1)  —  (^i+l  >  Qto;+2  )  •  •  •  >  )  (4.16) 

for  1  <  i  <  N .  An  iterative  approximation  algorithm  similar  to  the  one  developed  in  Section  4  is 
also  available  for  this  case;  details  are  omitted  for  the  sake  of  brevity. 

5.  NUMERICAL  EXAMPLES 

The  accuracy  of  the  algorithm  is  tested  against  simulations,  and  both  relative  and  absolute 
errors  in  the  approximations  are  listed.  Although  relative  errors  are  usually  a  more  important 
measure  of  accuracy  of  an  approximation  scheme,  absolute  errors  play  a  more  important  role 
when  determining  the  relative  or  absolute  errors  associated  with  the  computation  of  some  of  the 
performance  measures  of  interest.  To  illustrate  this  point,  average  queue  sizes  at  each  buffer  location 
are  calculated.  Both  the  relative  and  the  absolute  errors  for  this  performance  measure  are  better 
predicted  by  the  absolute  errors  in  the  calculation  of  the  probabilities.  For  instance,  in  Example 
1,  the  approximation  of  the  probability  of  having  an  empty  buffer  at  the  first  queue  yields  17.6% 
relative  error,  while  the  absolute  error  is  only  about  1%.  The  relative  error  in  the  average  queue 
size  for  this  buffer  is  1.9%.  Generally  speaking,  in  all  the  examples  the  approximations  are  within 
±0.02  of  the  simulations.  Such  an  accuracy  is  considered  reasonable  since,  due  to  blocking,  the 
system  usually  reaches  steady  state  rather  slowly.  Therefore,  the  simulations  are  also  believed  to 
be  accurate  only  up  to  the  second  digit  after  the  decimal  point. 

The  first  five  examples  are  for  the  discrete-time  model  under  the  immediate  blocking  strategy 
and  involves  only  geometric  servers.  Examples  6-17  are  for  the  continuous-time  model  operating 
under  the  non-immediate  blocking  policy  and  involves  more  general  service  time  distributions. 
The  accuracy  of  the  approximations  seems  not  to  be  affected  by  high  blocking  probabilities  (e.g., 
Examples  1,6,7,9-12)  or  in  the  presence  of  bottleneck(s)  (e.g.,  Examples  11  and  12).  Also,  the 
approximations  work  equally  well  even  when  service  time  distributions  have  small  or  large  squared 
coefficient  of  variations  (e.g.,  c2  =  0.2,  i  =  2,3  in  Example  10  and  2  <  c2  <  10  in  Examples  13-17). 

Examples  2-5  indicate  that  the  reversibility  property  discussed  in  Chapter  3  of  [6]  (see  also 
[3])  also  hold  for  a  general  tandem  line  under  the  immediate  blocking  strategy.  In  Examples  2  and 
3,  where  the  tandem  line  is  totally  symmetric,  reversing  the  order  of  the  servers  yield  the  exact 
same  probabilities.  Example  5  is  the  same  as  Example  4  except  that  the  order  of  both  the  servers 
and  the  intermediate  buffers  is  now  reversed.  The  approximation  scheme  yields  totally  symmetric 
probabilities.  In  the  simulation  results  for  these  examples,  the  corresponding  probabilities  change 
in  the  third  digit  after  the  decimal  point,  thereby  providing  an  indication  of  the  margin  of  accuracy 
of  the  simulations. 

For  the  numerical  values  of  the  Examples  8-13,  three  independent  simulation  runs  are  per¬ 
formed  on  the  Performance  Analysis  Workstation  (PAW)  [12].  The  results  obtained  from  each 
simulation  were  within  ±0.01  of  each  other.  The  results  listed  in  the  simulation  column  in  Tables 
5.8-5.12  are  the  arithmetic  average  of  these  three  simulation  runs.  In  these  examples,  Si  denotes 
the  service  time  distribution  of  server  i,  while  Er(r;  (.1)  (resp.  Er(2\  //1,  /t2))  denotes  the  r-stage 
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Erlang  (resp.  2-stage  Generalized  Erlang)  distribution  with  parameter(s)  pi  (resp.  /zi  and  H2),  i. 
e.,  is  the  average  time  spent  at  each  stage. 

Examples  14-17  are  given  in  order  to  compare  the  proposed  algorithm  (indicated  by  the  column 
GM)  against  the  one  proposed  by  Jun  and  Perros  (JP)  [10].  In  these  examples,  the  first  server 
generates  Poisson  arrivals  with  rate  A  and  jobs  that  arrive  into  the  system  when  the  first  buffer 
is  full  are  assumed  lost.  The  service  times  at  each  node  are  two-stage  Coxian  with  mean  l//i  and 
squared  coefficient  of  variation  c2.  The  corresponding  PH-representation  is  given  by 

^x2)  and  «  =  (1»°)’ 

where 

7i  =  2/i  ,  72  =  £  and  9  =  ^  ■ 

The  simulation  data  is  taken  from  [10],  where  the  average  queue  length  and  the  probabilities  of 
having  an  empty  buffer  and  a  full  buffer,  denoted  respectively  by  E(L),  p( 0)  and  p(N),  are  also 
available  for  each  buffer.  The  comparison  of  the  corresponding  data  shows  that  both  approximation 
algorithms  perform  well  relatively  to  the  simulation  data  and  have  comparable  accuracy  except  in 
Example  17  where  all  the  servers  have  high  variability.  In  this  example  the  algorithm  presented  in 
this  paper  provides  a  better  approximation  in  all  the  entries  of  Table  5.17.  However,  the  algorithm 
of  Jun  and  Perros  was  2-3  times  faster  than  the  proposed  algorithm.  This  is  mainly  because 
the  proposed  algorithm  is  implemented  for  general  PH-distributions  and  does  not  make  use  of  the 
special  structure  of  the  matrices  in  the  effective  representations.  As  mentioned  at  the  end  of  Section 
3,  while  several  steps  can  be  taken  to  reduce  the  CPU  time,  the  emphasis  in  the  present  paper 
is  mainly  on  the  structure  of  the  approximation  methodology  rather  than  on  the  implementation 
details. 
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Example  1  :  N  =  3,  K\  =  2,  if 2  =  1.  Immediate  blocking. 
Geometric  servers;  p,-  =  0.5,  i  =  1,2,3. 


Buffer 

Queue  Size 

Approx. 

Exact  Sohn.  Rel.  Error 

Abs.  Error 

1 

0 

0.0526 

0.0638 

0.1755 

0.0112 

1 

0.4211 

0.4256 

0.0106 

0.0045 

2 

0.5263 

0.5106 

-0.0307 

-0.0157 

Average  queue  length: 

1.4737 

1.4468 

-0.0186 

-0.0269 

2 

0 

0.5263 

0.5106 

-0.0307 

-0.0157 

1 

0.4737 

0.4894 

0.0321 

0.0157 

Average  queue  length: 

0.4737 

0.4894 

0.0321 

0.0157 

Table  5.1. 

Example  2  : 

N  -  4,  Ki  =  2,  i  =  1,2,3. 
Geometric  servers;  pi  =  0.5, 

Immediate  blocking. 
i  =  1,2, 3, 4. 

Buffer 

Queue  Size 

Approx. 

Exact  Sol’n. 

Rel.  Error 

Abs.  Error 

1 

0 

0.1444 

0.1652 

0.1259 

0.0208 

1 

0.4470 

0.4596 

0.0274 

0.0126 

2 

0.4086 

0.3754 

-0.0884 

-0.0332 

Average  queue  length: 

1.2643 

1.2104 

-0.0444 

-0.0538 

2 

0 

0.2643 

0.2621 

-0.0084 

-0.0022 

1 

0.4715 

0.4757 

0.0088 

0.0042 

2 

0.2643 

0.2621 

-0.0084 

-  -0.0022 

Average  queue  length: 

1.0001 

0.9999 

-0.0002 

-0.0002 

3 

0 

0.4086 

0.3754 

-0.0884 

-0.0332 

1 

0.4470 

0.4596 

0.0274 

0.0126 

2 

0.1444 

0.1652 

0.1259 

0.0208 

Average  queue  length: 

0.7358 

0.7900 

0.0686 

0.0542 

Table  5.2. 
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Example  3  :  N  —  4,  Ki  =  4,  t  =  1,2,3.  Immediate  blocking. 
Geometric  servers;  Pi  =  0.5,  i  =  1,2, 3, 4. 


Buffer  Queue  Size 

Approx. 

Simulation 

Rel.  Error 

Abs.  Error 

1  0 

0.074 

0.079 

0.063 

0.005 

1 

0.174 

0.186 

0.065 

0.012 

2 

0.234 

0.235 

0.004 

0.001 

3 

0.313 

0.305 

-0.026 

-0.008 

4 

0.205 

0.195 

-0.054 

-0.010 

Average  queue  length: 

2.401 

2.311 

-0.039 

-0.090 

2  0 

0.131 

0.127 

-0.031 

-0.004 

1 

0.246 

0.247 

0.004 

0.001 

2 

0.247 

0.247 

0.000 

0.000 

3 

0.246 

0.249 

0.012 

0.003 

4 

0.131 

0.130 

0.008 

0.001 

Average  queue  length: 

2.002 

2.008 

0.003 

0.006 

3  0 

0.205 

0.195 

-0.051 

-0.010 

1 

0.313 

0.308 

-0.016 

-0.005 

2 

0.234 

0.238 

0.017 

0.004 

3 

0.174 

0.182 

0.044 

0.008 

4 

0.074 

0.076 

0.026 

0.002 

Average  queue  length: 

1.599 

1.634 

0.021 

0.035 

Table  5.3. 
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Example  4  :  N  =  5,  K\  =  K3  —  2,  K2  =  A4  =  3.  Immediate  blocking. 

Geometric  servers;  pi  =  0.7,  p2  =  0.8,  p3  =  0.9,  p4  =  0.75,  p5  =  0.6. 


Buffer 

Queue  Size 

Approx. 

Simulation 

Rel.  Error 

Abs.  Error 

1 

0 

0.168 

0.176 

0.045 

0.008 

1 

0.599 

0.601 

0.003 

0.002 

2 

0.233 

0.223 

-0.045 

-0.010 

Average 

queue  length: 

1.065 

1.047 

-0.017 

-0.018 

2 

0 

0.109 

0.118 

0.076 

0.009 

1 

0.340 

0.334 

-0.018 

-0.006 

2 

0.389 

0.377 

0.032 

0.012 

3 

0.162 

0.171 

0.053 

0.009 

Average 

queue  length: 

1.604 

1.601 

-0.002 

-0.003 

3 

0 

0.116 

0.103 

-0.126 

-0.013 

1 

0.590 

0.589 

-0.002 

-0.001 

2 

0.294 

0.308 

0.045 

0.014 

Average 

queue  length: 

1.178 

1.205 

0.022 

0.027 

4 

0 

0.106 

0.091 

-0.165 

-0.015 

1 

0.331 

0.305 

-0.085 

-0.026 

2 

0.395 

0.416 

0.050 

0.021 

3 

0.168 

0.188 

0.106 

0.020 

Average 

queue  length: 

1.625 

1.701 

0.045 

0.076 

Table  5.4. 
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Example  5  :  N  =  5, 

II 

II 

CO 

II 

.  « 

Ki  =  2.  Immediate  blocking. 

Geometric 

servers;  p\  =  0.6,  p2 

=  0.75,  p3  =  0.9,  pA 

=  0.8,  p5  =  0.7. 

Buffer  Queue  Size 

Approx. 

Simulation 

Rel.  Error 

Abs.  Error 

1  0 

0.168 

0.191 

0.120 

0.023 

1 

0.395 

0.413 

0.044 

0.018 

2 

0.331 

0.303 

-0.092 

-0.028 

3 

0.106 

0.093 

-0.140 

-0.013 

Average  queue  length: 

1.375 

1.298 

-0.059 

-0.077 

2  0 

0.294 

0.308 

0.045 

0.014 

1 

0.590 

0.589 

-0.002 

-0.001 

2 

0.116 

0.103 

-0.126 

-0.013 

Average  queue  length: 

0.822 

0.776 

-0.059 

-0.046 

3  0 

0.175 

0.172 

0.058 

0.010 

1 

0.389 

0.376 

-0.035 

-0.013 

2 

0.340 

0.333 

-0.021 

-0.007 

3 

0.109 

0.119 

0.084 

0.010 

Average  queue  length: 

1.396 

1.337 

-0.044 

-0.059 

4  0 

0.230 

0.224 

-0.040 

-0.009 

1 

0.559 

0.604 

0.008 

0.005 

2 

0.168 

0.172 

0.023 

0.004 

Average  queue  length: 

0.935 

0.939 

0.004 

0.004 

Table  5.5. 
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Example  6  :  N  =  3,  K{  =  1,  i  =  1,2.  Non-immediate  blocking. 


Exponential 

servers;  /x,  = 

1,  i  =  1,2,3. 

Buffer 

Queue  Size 

Approx. 

Exact  Sol’n. 

Rel.  Error 

Abs.  Error 

1 

0 

0.222 

0.205 

-0.083 

-0.017 

1 

0.778 

0.795 

0.021 

0.017 

2 

0 

0.444 

0.436 

-0.018 

-0.008 

1 

0.556 

0.564 

0.014 

0.008 

Table  5.6. 

Example  7  : 

N  =  3,  Ki  = 
Exponential 

1,  *  =  1,2. 

servers;  mi  - 

Non-immediate  blocking. 

=  1,  M2  =  1-75,  M3  =  1.5 . 

Buffer 

Queue 

Size 

Approx. 

Exact  Sol’n. 

Rel.  Error 

Abs.  Error 

1 

0 

0.401 

0.398 

1 

0.599 

0.602 

2 

0 

0.509 

0.519 

0.010 

1 

0.491 

0.481 

-0.010 

Table  5.7. 

Example  8 

:  N  —  3,  K\  =  2,  I{2  =  3.  Non-immediate  blocking. 

Si  =  Er{ 2;  2, 1),  S3  =  Er( 2;  2.5, 1), 

5*2  =Hyperexponential;  Q  =  diag(— 1,  — 0.5),  a  =  (0.3, 0.7). 

Buffer 

Queue  Size  Approx. 

Simulation 

Rel.  Error 

Abs.  Error 

1 

0  0.158 

0.155 

-0.019 

-0.003 

1  0.242 

0.248 

0.024 

0.006 

2  0.601 

0.597 

-0.007 

-0.004 

Average  queue  length:  1.444 

1.442 

-0.001 

-0.002 

2 

0  0.355 

0.351 

-0.011 

-0.004 

1  0.286 

0.288 

0.007 

0.002 

2  0.185 

0.195 

0.051 

0.010 

3  0.174 

0.166 

-0.048 

-0.008 

Average  queue  length:  1.178 

1.176 

-0.002 

-0.002 

Table  5.8. 
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Example  9  :  N  =  3,  K\  =  3,  K2  =  2.  Nonimmediate  blocking. 

Si  =  Er( 2;  1),  S2  =  Er(2;  &),  S3  =  £r(2;  1,  &). 


Buffer 

Queue  Size 

Approx. 

Simulation 

Rel.  Error 

Abs.  Error 

1 

0 

0.035 

0.027 

-0.296 

-0.008 

1 

0.104 

0.091 

-0.143 

-0.013 

2 

0.205 

0.207 

0.010 

0.002 

3 

0.656 

0.675 

0.028 

0.019 

Average  queue  length: 

2.482 

2.530 

0.019 

0.048 

2 

0 

0.235 

0.231 

-0.017 

-0.004 

1 

0.353 

0.353 

0.000 

0.000 

2 

0.412 

0.415 

0.007 

0.003 

Average  queue  length: 

1.177 

1.183 

0.005 

0.006 

Table  5.9. 

Example  10  : 

:  N  =  3,  K\  —  2,  K2  —  3.  Non-immediate.blocking. 

S\  =  Er(  1;  0.25),  S2  =  Er( 5;  1.5),  S3  =  Er( 5;  1). 

Buffer 

Queue  Size 

Approx. 

Simulation  Rel.  Error 

Abs.  Error 

1 

0 

0.171 

0.170 

-0.006 

-0.001 

1 

0.269 

0.251 

-0.072 

-0.018 

2 

0.560 

0.579 

0.033 

0.019 

Average  queue  length: 

1.389 

1.397 

0.006 

0.008 

2 

0 

0.042 

0.058 

0.276 

'  0.016 

1 

0.126 

0.129 

0.023 

0.003 

2 

0.240 

0.225 

-0.067 

-0.015 

3 

0.592 

0.588 

-0.007 

-0.004 

Average  queue  length: 

2.382 

2.343 

-0.017 

-0.039 

Table  5.10. 
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Example  11  :  N  =  6,  Ki  =  3,  1  <  i  <  5.  Non-immediate  blocking. 

Si  =  Er( 2;  1),  1  <  *  f  4  <  6,  S4*  =  Er{ 2;  0.5). 


Average  Queue  Lengths 


Buffer _ Approx. _ Simulation _ Rel.  Error 


1 

2.66 

2.85 

0.067 

2 

2.79 

2.84 

0.018 

3 

2.82 

2.83 

0.004 

4 

0.68 

0.67 

-0.015 

5 

0.68 

0.69 

-0.014 

Table  5.11. 

Example  12  :  N 

=  10,  Non-immediate  blocking. 

<  Exponential,  /i  =  0.5, 

i  =  1, 3, 6, 8 

f 

3,  i  —  1,2, 5, 6, 8 

Exponential,  /x  =  1, 

i  =  10 

II 

O' 

*■*< 

2,  *  =  4,7,9 

Si  = 

Exponential,  /x  =  0.2, 

i  =  4* 

l 

<?■». 

II 

w 

Er{  2;  1), 

i  =  2, 5,  7 

l  Er(2;  0.5), 

i  =  9*. 

Average  Queue  Lengths 

Buffer  Approx. 

Simulation  Rel.  Error 

1 

2.69 

2.84  0.053 

2 

2.80 

2.87  0.024 

3 

3.77 

3.77  0.000 

4 

0.52 

0.52  0.000 

5 

0.58 

0.63  0.063 

6 

0.66 

0.77  0.143 

7 

0.72 

0.77  0.065 

8 

1.59 

1.58  -0.006 

9 

0.21 

0.21  0.000 

Table  5.12. 


*  Bottleneck 
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Example  13  :  N  =  3,  K\  =  K?  =  3.  Non-immediate  blocking. 

<2,-  =  diag(-0.1, -10),  a;  =  (0.5, 0.5),  c?  =  2.93,  i—  1,2,3  . 
(Hyperexponential) 


Buffer 

Queue  Size 

Approx. 

Simulation 

Rel.  Error 

Abs.  Error 

1 

0 

0.214 

0.205 

-0.044 

-0.009 

1 

0.104 

0.114 

0.088 

0.010 

2 

0.118 

0.122 

0.033 

0.004 

3 

0.563 

0.559 

-0.007 

-0.004 

Average  queue  length: 

2.029 

2.019 

-0.005 

-0.010 

2 

0 

0.428 

0.427 

-0.002 

-0.001 

1 

0.135 

0.133 

-0.015 

-0.002 

2 

0.118 

0.118 

0.000 

0.000 

3 

0.319 

0.322 

0.009 

0.003 

Average  queue  length: 

1.328 

1.335 

0.005 

0.007 

Table  5.13. 

Example  14 

:  N  =  3,  I<i 

II 

CO 

II 

II 

3.  Non-immediate  blocking. 

A  =  3.0,  pLi 

—  4,  *  —  1,2, 3,  cf 

—  2,  c\  —  4,  c|  — 

8. 

Buffer 

Simulation 

GM 

Rel.  Error  JP 

Rel.  Error 

1 

p(  o) 

0.253 

0.264 

-0.043 

0.244 

0.036 

p(K) 

0.353 

0.343 

0.028 

0.361 

-0.023 

E(L ) 

1.646 

1.609 

0.022 

1.668 

.  -0.013 

2 

p(°) 

0.336 

0.345 

-0.027 

0.337 

-0.003 

p(K) 

0.374 

0.364 

0.027 

0.386 

-0.032 

E{L) 

1.537 

1.507 

0.020 

1.559 

-0.014 

3 

v{  0) 

0.516 

0.518 

-0.004 

0.521 

-0.010 

p(K) 

0.239 

0.236 

0.013 

0.250 

-0.046 

E(L) 

1.053 

1.045 

0.008 

1.069 

-0.015 

Table  5.14. 
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Example  15  :  N  =  3,  K\  =  I(2  =  ^3  =  5.  Non-immediate  blocking. 


II 

0 

Pi  =  4,  c?  =  2, 

i  =  1,2,3. 

Buffer 

Simulation 

GM 

Rel.  Error 

JP 

Rel.  Error 

1 

K°) 

0.110 

0.118 

-0.073 

0.116 

-0.055 

P(K) 

0.310 

0.303 

0.023 

0.310 

0.000 

E(L) 

3.125 

3.094 

0.010 

3.088 

0.012 

2 

p(0) 

0.193 

0.199 

-0.031 

0.197 

-0.021 

P(K ) 

0.325 

0.323 

0.006 

0.319 

0.018 

E(L) 

2.822 

2.796 

0.009 

2.798 

0.009 

3 

p(  o) 

0.314 

0.316 

-0.006 

0.310 

0.013 

p(K) 

0.201 

0.201 

0.000 

0.198 

0.015 

E(L) 

2.108 

2.101 

0.003 

2.115 

-0.003 

Table  5.15. 

Example  16 

:  N  =  3, 
A  =  4.0, 

Ki  =  K3  =  5,  I(2 

Hi  =  P3  =  3,  H2 

=  3. 

=  5, 

Non-immediate  blocking. 
c?  =  4,  *=1,2,3. 

Buffer 

Simulation 

GM 

Rel.  Error 

JP 

Rel.  Error 

1 

JP(°) 

0.059 

0.070 

-0.186 

-0.119 

p(K) 

0.495 

0.478 

0.034 

0.495 

0.000 

E(L) 

3.763 

3.333 

0.114 

3.711 

0.014 

2 

p(  0) 

0.378 

0.383 

-0.013 

0.386 

-  -0.021 

p(K) 

0.375 

0.363 

0.032 

0.367 

0.021 

E(L) 

1.486 

1.455 

0.021 

1.464 

0.015 

3 

p(  0) 

0.327 

0.326 

0.003 

0.326 

0.003 

p(K) 

0.298 

0.295 

0.010 

0.296 

0.007 

E(L) 

2.373 

2.369 

0.002 

2.388 

0.006 

Table  5.16. 
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Example  17  :  N  =  3,  K\  =  if 2  =  if 3  =  3.  Non-immediate  blocking. 
A  =  4.0,  /i»  =  4,  c?  =  8,  i  =  1,2,3. 


Buffer 

Simulation  GM 

Rel.  Error  JP 

Rel.  Error 

1 

p(  0) 

0.187 

0.196 

-0.048 

0.172 

0.080 

p(/0 

0.524 

0.514 

0.019 

0.540 

-0.031 

E(X) 

2.004 

1.969 

0.017 

2.051 

-0.023 

2 

p(o) 

0.333 

0.338 

-0.015 

0.340 

-0.021 

K/O 

0.438 

0.431 

0.016 

0.469 

-0.071 

E{L) 

1.647 

1.626 

0.013 

1.692 

-0.027 

3 

p(  0) 

0.532 

0.528 

0.008 

0.540 

-0.015 

p(K) 

0.254 

0.256 

-0.008 

0.272 

-0.071 

E(L ) 

1.060 

1.071 

-0.010 

1.088 

-0.026 

Table  5.17. 
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APPENDIX 


THE  TWO  NODE  TANDEM  SYSTEM 
The  discrete-time  model 

Consider  the  tandem  system  described  in  Section  2  when  N  =  2  and  K\  =  K.  A  natural  state 
space  E  for  this  system  is  the  one  that  contains  r  :=  ( K  —  l)mim2  +  m\  +  m2  states  with 


'  (i,  0)  ,  1  <  i  <  mi  ,  k  =  0, 

E  —  (i,  fc,  j) ,  1  <  i  <  mi,  0  <  k  <  K  and  1  <  j  <  m2, 
.  (K,j)  ,  1  <j<m2,  k-I(, 


where  k  indicates  the  buffer  size,  and  i  and  j  represent  the  service  phase  in  the  first  and  the  second 
node  server,  respectively. 

A  discrete-time  Markov  chain  is  naturally  associated  with  E;  its  one-step  transition  matrix  T 
is  described  in  [6].  Any  invariant  probability  vector  for  T  is  denoted  by  the  1  x  r  row  vector  n,  which 
is  partitioned  into  K  +  1  blocks  of  components,  say  x  =  (7r0 ,  tti,  . . . ,  xk),  with  x0,Xk,  0  <  k  <  K, 
and  xk  being  row  vectors  of  dimension  1  x  mi,  1  X  mi  m2  and  1  x  m2,  respectively. 

First,  several  matrices  are  defined  in  terms  of  the  system  parameters. 


-M  —  (fmi  <*1)  ®  fnij  d~  (^mj  ^1  Q l)  ®  Q 2 

A  =  An ,  ®  (An2  ^2)  4"  Ql  ®  (®m2  Q2) 

R  =  M  N~l 


S  =  [(Imi-Qi)®a2]  N-1 
U  =  {Ql®  Pl )  {Im2  ~  Ql)~l 


W  =  Imi  + 


(s 


K- 1 
k=l 


S  Rk' 


(-fmj  ®  ^7n2  ) 


Z  =  W  +  SRk~ 2  U  em2  x 


mi  m2  X  mi  m2, 
mi  m2  X  mi  m2, 

mi  m2  x  mi  m2, 

mi  X  mi  m2, 
mi  m2  X  m2, 

mi  x  mi, 
mi  x  mi, 


where  the  1  X  mi  row  vector  a;  =  «i(fmi  —  Qi)_1/ai(/m,  —  <5i)_1em,  is  the  invariant  probability 
vector  of  Qi  +  Pi ■ 

When  K  >  1,  if  either  the  matrix  M  or  N  is  invertible,  then  a  closed-form  solution  can  be 
obtained  for  the  vector  x.  The  main  result  of  [6]  is  summarized  in  the  following  theorem  when  N 
is  invertible. 


Theorem  A.l.  If  the  matrix  N  is  invertible,  then  any  invariant  probability  vector  x  —  (7To, . . . ,  Xj<) 
has  the  following  matrix  geometric  form 

xk  —  xq  S  Rk~l  ,  1  <  k  <  K  (A. la) 

and 

xK  =  t t0  S  RK~ 2  U  ,  (A. 16) 


where  the  vector  xq  satisfies  the  equation 


When  K  =  1,  i.  e.,  there  is  no  intermediate  buffer,  a  simpler  expression  is  available  in  [4,  6]. 

The  following  necessary  and  sufficient  conditions  for  the  invertibility  of  the  matrices  M  and  N 
are  given  in  [6].  To  state  them,  let  C  be  the  open  disc  centered  at  (|,0)  in  the  complex  plane  with 
radius  \  and  let  Sp(X)  denote  the  spectrum  of  the  matrix  X. 

Lemma  A.l.  If  Sp(A)  C  C  (resp.  Sp(B )  C  C)  then  the  matrix  N  (resp.  M)  is  nonsingular. 
Lemma  A. 2.  The  matrix  N  (resp.  M)  is  singular  if  the  matrix  A  (resp.  B)  is  singular. 


Let  u  (resp.  v)  be  the  steady  state  probability  vector  of  the  Markov  chain  associated  with  E 
embedded  at  points  of  arrival  (resp.  departure).  These  vectors  are  partitioned  into  K  blocks  of 
components,  say  u  =  (uo,  «i, . . . ,  ujc_i)  and  v  =  (u0,  V\, . . . ,  uj<-_i),  with  uk  and  vk ,  0  <  k  <  K, 
being  1  X  m2  and  1  X  mi  row  vectors,  respectively.  The  Ith  component  of  uk  (resp.  vk)  is  the  joint 
steady-state  probability  that  at  a  service  completion  epoch  at  the  first  (resp.  second)  server,  there 
are  k  jobs  in  the  buffer,  excluding  the  job  that  has  just  been  serviced,  and  that  the  phase  of  the 
second  (resp.  first)  server  is  i.  The  following  theorem  establishes  simple  relations  between  u,  v  and 
7r;  the  proofs  are  similar  to  the  ones  given  in  [3]  and  are  therefore  omitted. 

Theorem  A. 2.  Define  the  scalars  c  and  d  by 


and 


Ii- 1 

c  =  7T0pi  +  Y  rk(Pi  ®  e  m2  ) 

1 

K- 1 

d  =  ni<P2  +  Y  7rfc(e  771 1  ®Pi)  ■ 
l 


(A.3a) 


(A. 36) 


For  K  >  1  the  relations 

cuq  =  (7ToPi)a2  +  7Ti(pi  ®p2a2) ,  (A.4a) 

cuk  =  7Tfc(pi  <g>  Q2)  +  7Tfc+1(pi  ®p2a2),  1  <k  <  K  -  1,  (A.46) 

CUk-I  =  7T/f_i(pi  ®Q2 )  (A.4c) 

and 

dv0  =  wi(Qi  ®  p2) ,  (A. 5a) 

dvk  =  tt/c(piOi  ®p2)  +  nk+1(Qi  ®  p2)  1  <  k  <  K  -  1,  (A. 56) 

dvK- 1  =  tt/v  — ! (Pi Qi  ®  p2)  +  (*■ kP2)<*i  ,  (A. 5c) 


hold  true,  while  for  I(  =  1,  Uq  =  a2  and  v0  =  «i- 
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The  continuous-time  model 

Equations  (A.l)  and  (A.2)  of  Theorem  A.l  also  holds  for  the  continuous-time  formulation  of 
the  two  node  model.  However,  the  matrices  M,  N,  S  and  U  are  now  given  by 

M  =  ( Imi  -  emi  ax)  ®  Q2  +  Qi  <8  Im2  mi  m2  x  mi  m2, 

N  =  Qi®  (/m2  -  em2  a2)  +  Imi  ®  Q2  mi  m2  X  mx  m2, 

S  =  [Q i  ®  a2]  A-1  mi  X  mi  m2, 

A  =  —  (pi  ®  Im2)Q2l  mi  m2  X  m2 

and  the  matrices  M  and  N  are  always  both  invertible. 

For  the  continuous-time  model,  the  equations  (A.4)  and  (A. 5)  of  Theorem  A.2  become 


cu0  =  (ir0pi)a2  ,  {A. 6a) 

cuk  -  nk(pi  ®  /m2),  1  <  k  <  K,  (A.6b) 


and 


dvk  =  TTfc+i (/mi  ®p2)  0  <  &  <  K  -  1, 


(A. 7a) 


du/c_i  =  (7T/<-p2)ai  , 


(A.76) 


where  the  scalars  c  and  d  are  again  given  by  (A.3). 
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